Chebyshev-BdG: an efficient numerical approach to inhomogeneous superconductivity 
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We propose a highly efficient numerical method to describe inhomogeneous superconductivity by 
using the kernel polynomial method in order to calculate the Green's functions of a superconductor. 
Broken translational invariance of any type (impurities, surfaces or magnetic fields) can be easily 
incorporated. We show that limitations due to system size can be easily circumvented and therefore 
this method opens the way for the study of scenarios and/or geometries that were unaccessible 
before. The proposed method is highly efficient and amenable to large scale parallel computation. 
Although we only use it in the context of superconductivity, it is applicable to other inhomogeneous 
mean-field theories. 
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In the past decades the mean-field description of inho- 
mogeneous superconductivity through the Bogoliubov- 
de Gennes (BdG) equations has been highly successful 
in uncovering novel phenomena. Because in the pres- 
ence of broken translational invariance one needs to use 
a real space formulation, the numerical simulation be- 
comes computationally involved. While alternative ap- 
proaches to inhomogeneous superconductivity like qua- 
siclassical approximations or Ginzburg-Landau methods 
exist, the need for a fully quantum mechanical approach 
has become imperative. This is manifest in questions re- 
garding high-Tc superconductors for which the supercon- 
ducting coherence length is of the order of the Fermi wave 
length, or in questions regarding nanoscale superconduc- 
tivity for which the superconducting coherence length is 
comparable to the system size. 

The BdG equations have been extensively used in a 
multitude of situations where translational symmetry is 
broken. Examples include the description of quasipar- 
ticles in s-wave or d-wave vortices |lH3|, self-consistent 
calculation of order parameters (OP) and local density 
of states (LDOS) near surfaces and interfaces y-|9|, self- 
consistent description of magnetic and non-magnetic im- 
purities in superconductors 10l4l2J. calculation of DC 
Josephson currents through weak links [5|, |6| , uncovering 
of the effect of electron confinement on superconductivity 
[H, etc. 

Throughout these studies several methods of solving 
the BdG equations have been employed. First, after 
a discretization of the mean-field Hamiltonian one can 
use the straightforward approach of diagonalizing exactly 
the resulting Hamiltonian. Although exact diagonaliza- 
tion can in principle treat any inhomogeneous situation 
it has severe limitations on the size of the discretization 
grid. One cure is to recover translational symmetry ei- 
ther by considering surfaces and interfaces or by consider- 
ing highly symmetric geometries (cylindrical or square). 
This way, by using a Fourier transformation in the direc- 
tion which retains the translational invariance, one can 
reduce the dimensionality of the problem: for each value 



of the momentum vector, we have to solve the BdG equa- 
tions of dimension d — 1, where d is the dimension of 
the initial problem. Another way of circumventing the 
size limitations of the exact diagonalization is the use 
of super-cells. This is done by considering an inhomo- 
geneous finite size region which is then replicated in all 
directions. The super-cell method is thus able to de- 
crease the spacing between eigen- energies and obtain a 
much smoother LDOS. This is again achieved by diago- 
nalizing the Hamiltonian of the finite size region for each 
momentum vector defined by the super-cell lattice. 

A completely different approach is based on approx- 
imating the Green's functions. In this case the eigen- 
energies will appear as poles of the Green's function while 
the wave-functions amplitudes will appear as weights of 
the poles. One such method is the recursive method 
based on the Lanczos procedure [9|, llj| . The approach we 
use here is similar in spirit but has several benefits when 
compared to the recursive method. We will show how 
the Green's function can be efficiently expanded in series 
of Chebyshev polynomials. The paper is organized as 
follows: first we will introduce a general model Hamilto- 
nian which is typically used for describing inhomogeneous 
superconductors. We will next present the Chebyshev- 
Bogoliubov-de Gennes (CBdG) method and show, by an 
example, how this method can be implemented. 

The Bogoliubov-de Gennes equations are mean-field 
coupled equations which describe the behavior of elec- 
trons and holes in superconductors. If we consider second 
quantization and work within the Nambu spinor formal- 
ism, a general Hamiltonian describing superconductivity 
can be written as follows: 
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where T-Lij is a 2 x 2 matrix: 
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Ci describes an on-site potential due to impurities, ^ is the 
chemical potential, Uj describes hopping between near- 
est neighbor sites while Aj(Aij) are the on-site(nearest 
neighbor) superconducting order parameters. The ef- 
fect of a magnetic field is contained in the complex or- 
der parameters through the usual Peierls phases i^ = 
\tij\ exp(i -^ J^ Aijdl), where Aij is the vector potential 
and 4>Q = h/2e is the flux quantum. 

The quantity of interest is the 2x2 Green's function, 
which is defined as: 
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its spectrum is contained in the [—1,1] interval. We 
therefore have to work with the rescaled Hamiltonian 
Ti = {% — 16)/a and rescaled energies E = [E — b)/a, 
cj = (w - b)/a where a = {Emax ~ £^rmn)/(2 - t]) and 
b ~ {Ejnax + Emin)/'^, whcrc Tj > Q IS. & small number. It 
is not essential to have accurate bounds on the spectrum, 
thus a quick Lanczos procedure to find Emax and Emm 
can be used. 

If we consider the regular Green's function in the 
Lehman representation we can write for its imaginary 
part: 

-sGljiuj + zrj) = -7:J2{cn\k){k\c]^)SiuJ - E^) (9) 



where G{u!+i7]) = [uj+ii]—'H]^^ and \vac) is the vacuum . 
The diagonal and off-diagonal components are the normal 
and anomalous Green's functions: 
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where |cL) = cL|wac) creates a spin-up electron and 
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\cii) = Cii\vac) destroys a spin-down electron. For fi- 
nite temperatures the expectation value also contains a 
thermal average. 

As mentioned before, the Green's function can be ap- 
proximated by using a Lanczos procedure to invert the 
Hamiltonian [9|, [ij] . This method has proven to be effi- 
cient mostly in the homogeneous case or when the Lanc- 
zos procedure can be easily extrapolated. The need of 
extrapolation is of utmost importance because due to 
numerical round-off errors the Lanczos procedure is un- 
stable after a number of iterations. Re-orthogonalization 
schemes exist but the method becomes less and less effi- 
cient. 

We therefore propose another approach to approxi- 
mate the Green's function. Our method is based on the 
Kernel Polynomial Method [15[ which expands the sin- 
gle particle Green's function into a series of Chebyshev 
polynomials. Any integrable function f{x) : [— 1, 1] ^ K, 
can be expanded as: 
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where Tn{x) = cos[n arccos(a;)] are the Chebyshev poly- 
nomials of first kind and (5o,n is the Kronecker delta func- 
tion. They are described by the following recursive rela- 
tions 



Tn+i{x) = 2xTn{x) - T„_i(x). 



(8) 



In order to be able to expand the Green's function, 
one needs first to rescale the Hamiltonian such that 



where {\k)} arc the eigenvectors and E'j. arc the eigenval- 
ues. If we now use the Chebyshev expansion we write: 
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where the coefficients can be calculated as the matrix 
elements of the Chebyshev polynomial of order n of "H : 
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With the use of a Kramcrs-Kronig relation, the real part 
of the Green's function can now be expanded in terms of 
Chebyshev polynomials of second kind Un{x) = sin[(n + 
1) arccos(a;)]/sin[arccos(a;)] [15| : 

oo 

5RGi/(cD + z,;) = -2 ^ a,V(*,.7)C^n(c^) (12) 

n=l 

After combining Eq. pO)) and Eq. (IT2t , the fuU Green's 
function can be written as: 
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withamij) = (c»t|T„(H)|c]^)/(l+Jo,n)- The procedure 
of finding the anomalous Green's function is identical, 
only the coefficients will be modified accordingly: 

««'(*. j) = (cL|T„(H)|ct )/(! + So,n) (14) 



The most important part of the calculation has now 
shifted to the calculation of the expansion coefficients 
a"^(i, j). Fortunately, due to the recurrence relation be- 
tween Chebyshev polynomials, see Eq. (|S]), these mo- 
ments can be obtained efficiently through a recursive pro- 
cedure. 

If wc define |j„) = r„('H)|cJ ), then after using the 
recursive property of Chebyshev polynomials [8] we can 
write: 
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At each iteration 



where |jo) 

step a]^{i,j) = (a|j„), where (1| = {c,^\ and (2| = (cj^|. 
It is important to note at this point that in the recur- 
sion defined by Eq. |15] the most intensive computation 
is a sparse matrix - vector muhiphcation. Moreover, the 
Hamiltonian matrix does not have to be stored since it al- 
ways has the same form, thus aUowing for simple rules for 
the multiplication. Another great benefit of this method 
is the possibility of obtaining in a single iteration all the 
normal and anomalous Green's functions, G\°'{u]), for all 

{i} and {a} when the starting vector is |cL)- As ex- 
plained in Ref. Il5l because we can only keep a finite 
number of terms in the expansion, one needs to convo- 
lute the approximated function with kernel polynomials 
in order to remedy the effect of Gibbs oscillations. This 
is imperative when approximating Green's functions be- 
cause of their discontinuous nature; the imaginary part 
is a summation over delta functions. We will use the 
Lorcntz kernel [l5|, since it allows for the manipulation 
of a Lorentzian broadened delta function. The expan- 
sion has the same form, but the coefficients have to be 
multiplied by factors defined by the Lorentz kernel: 
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where N is the total number of terms in the expan- 
sion and A is a real number. If we write the Lorentzian 
approximation as 5{x) = l/Trlinij^o ^/(s;^ + e^), then 
there is a direct relation between the broadening e and 
A: e = A/TV. This allows for a good control over the 
broadening of the Green's function's features, whether 
used artificially at zero temperature or naturally at fi- 
nite temperature. As we will show later, in certain situa- 
tions where interference between parts of the considered 
system is important, one needs a large number of coeffi- 
cients in order to accurately obtain the Green's function. 
In that case the only way to keep the broadening constant 
is by changing A accordingly. 

Once the Green's functions arc known, it is straight- 
forward to calculated physically relevant quantities. The 
local density of states can be calculated as: 



1 



-5a 



11(22) 



{E). 



The electron density is: 

[N'^iE, i) + N^{E, i)] f{E)dE. 



(17) 
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The order parameter, A^ — Uij (cifC^j,) is: 

A„- = ^C/y / ' Glf{E){l - 2f{E))dE, (19) 

where Ec is a cutoff energy (Debye energy for conven- 
tional superconductors or the bandwidth for cuprates). 
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FIG. 1. LDOS and a^ at the surface of a planar s-wave 
superconductor/normal metal system. L^ja = 300, L^ /a = 
20 and Ly/a = 500. 



The current density between grid points i and j is: 
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One of the great benefits of this method is that the 
Green's function is calculated separately for each grid 
point thus allowing for a trivial parallel implementation. 
An iteration can be started on a separate CPU for each 
grid point with a given order parameter profile. Next 
the order parameter for that grid point is updated in the 
Hamiltonian in order to achieve self-consistency. The 
method is general and it can be applied not only to 
any mean-field Hamiltonian but also to more complex 
band structures, multi-band superconductivity and even 
to three dimensional systems. Of course the number of 
operations increases dramatically but the calculation can 
be done even on a desktop computer since the Hamilto- 
nian is sparse. 



As an example we will show how the LDOS depends on 
the number of Chebyshev coefficients. We consider first 
a planar system composed of a normal metal of length 
L^ = 380a and an s-wave superconductor of length 
Ll = 20a while Ly = 500a. In Fig. [^a) we plot the 
LDOS at the surface of the normal metal region. Choos- 
ing A = O.li such that ^ « 7a, we observe in the LDOS 



s-wave/nortnal metal/s-wave normal metal 




(d) E/A=0.075 




(e)E/A=0.175 



X^ 



(f) E/A=0.425 



J L. 



- 0.1 



- 



-0.1 



50 60 70 80 90 10050 60 70 80 90 100 

X X 



FIG. 2. Local density of states around an impurity in a nor- 
mal metal for various energies. For panels (a)-(c) the normal 
metal is sandwiched between two s-wave superconducting re- 
gions located at a; < 50a and x > 100a respectively, while 
for panels (d)-(f) the whole system is normal. The LDOS for 
the corresponding clean systems (no impurity) is subtracted 
in both cases, for clarity. 



sion method based on the Lanczos method fails in these 
situations. 

To illustrate the power of the method we show in 
Fig. [2] the LDOS for a s-wave SC/normal metal/s- 
wave SC of size 50a/50a/50a x 100a in the presence 
of a non-magnetic impurity in the normal region Vi = 
yexp[— (ri — rio)^/a] with V = 2t and r^o = (60a, 50a). 
The left panels show the LDOS around the impurity for 
various sub-gap energies while the right panels show the 
LDOS for a homogeneous normal system. Modifications 
of the LDOS induced by the impurity are seen in both 
cases but for the s- wave/normal metal/s-wave system ex- 
tra states are induced by the interference of quasiparti- 
cles undergoing Andreev reflection at the superconduc- 
tor/normal metal interface and specular reflection at the 
impurity site. Andreev states of the clean multilayer sys- 
tem are destroyed by the impurity but new states appear 
due to impurity scattering. 

In conclusion we have introduced and demonstrated a 
new method of solving the mean- field self-consistent BdG 
equations by expanding the Nambu Green's functions in 
terms of Chcbyshcv polynomials. Because the method 
is stable the results are arbitrarily accurate since the ac- 
curacy is given by the number of moments kept in the 
expansion. The most expensive numerical operation is a 
sparse matrix - vector multiplication, thus allowing for 
large sized systems to be solved with little memory re- 
quirements. Moreover, since each grid point is calculated 
separately the method is amenable to trivial parallel im- 
plementations. The present method can be easily ex- 
panded to consider complex band structures, multi-band 
superconductivity, three dimensional system and other 
mean- field Hamiltonians. 
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the appearance of Andreev bound states below the super- 
conducting gap. In Fig. [ijb) we plot the moments of the 
Chebyshev expansion for three iteration sequences. Here 
we choose a constant broadening e ~ O.OOlt, thus the co- 
efficient A = t/N will modify the Chebyshev moments for 
each sequence. We observe an oscillatory behavior of the 
Chebyshev moments which is given by the interference of 
quasiparticles scattering off the normal/superconducting 
and normal/vacuum interfaces. Note that the Cheby- 
shev iteration is equivalent to a propagation of a quasi- 
particle defined by the starting vector |z). Interestingly 
the LDOS is not converged within e for A^ = 2000, in- 
stead a larger number of moments is needed. It is exactly 
for these type of systems that a stable method is essen- 
tial. When interference between quasiparticles scattered 
of distant regions of the system is important, an accurate 
solution requires a large number of moments. The recur- 
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